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A  MINIMUM-RESIDUAL  FINITE  ELEMENT  METHOD  FOR  THE 
CONVECTION-DIFFUSION  EQUATION 

JESSE  CHAN,  JOHN  A.  EVANS 

Abstract.  We  present  a  minimum- residual  finite  element  method  for  convection-diffusion  problems  in  a 
higher  order,  adaptive,  continuous  Galerkin  setting.  The  method  borrows  concepts  from  both  the  Discontinuous 
Petrov- Galerkin  (DPG)  method  by  Demkowicz  and  Gopalakrishnan  [1]  and  the  method  of  variational  stabilization 
by  Dahmen  et  al.  [2] ,  and  it  can  also  be  interpreted  as  a  variational  multiscale  method  in  which  the  fine-scales  are 
defined  through  a  dual-orthogonality  condition.  A  key  ingredient  in  the  method  is  the  proper  choice  of  norm  used 
to  measure  the  residual,  and  we  present  two  alternatives  which  are  observed  to  be  robust  in  both  convection  and 
diffusion-dominated  regimes.  Numerically  obtained  convergence  rates  are  given  in  2D,  and  benchmark  numerical 
examples  in  all  space  dimensions  are  shown  to  illustrate  the  behavior  of  the  method. 


1.  Introduction.  It  is  well  known  that  the  standard  Galerkin  finite  element  method  per¬ 
forms  very  poorly  for  the  convection-diffusion  equation  -  in  the  convection-dominated  case,  it 
experiences  spurious  oscillations  in  the  solution  that  grow  as  e  -A  0.  The  problem  is  connected 
back  to  a  loss  of  discrete  coercivity  with  respect  to  the  standard  H 1  norm  [3].  The  concept 
of  stabilized  methods  was  introduced  in  order  to  combat  such  oscillations,  the  most  popular  of 
which  is  the  Streamline  Upwind  Petrov- Galerkin  (SUPG)  method  [4].  The  method  can  be  in¬ 
terpreted  as  adding  a  sufficient  amount  of  artificial  viscosity  in  the  streamline  direction  in  order 
to  restore  discrete  coercivity  with  respect  to  a  new  “streamline-diffusion”  norm  [5].  SUPG  is 
also  an  example  of  a  residual-based  stabilization,  where  the  stabilization  mechanism  disappears 
as  the  strong  residual  of  the  equation  is  satisfied. 

A  connection  can  be  drawn  between  residual-based  stabilized  methods  and  Petrov- Galerkin 
schemes,  where  the  trial  (approximating)  functions  and  test  (weighting)  functions  are  allowed 
to  differ.  Specifically,  the  SUPG  method  can  be  interpreted  as  a  modification  of  standard  test 
functions1,  biasing  them  in  the  upwind  direction  based  on  mesh  size,  order  of  polynomial  approx¬ 
imation,  the  magnitude  of  convection,  and  the  diffusion  parameter.  More  recently,  the  concept 
of  Petrov- Galerkin  methods  has  been  connected  to  a  novel  minimum  residual  framework  -  it 
is  shown  that  the  minimization  of  a  specific  residual  corresponding  to  a  variational  formulation 
naturally  leads  to  the  concept  of  optimal  test  functions  [6]. 

Optimal  test  functions  resulting  from  residual  minimization  were  first  implemented  by 
Demkowicz  and  Gopalakrishnan  in  [1,  7].  The  connection  between  stabilization  and  least 
squares/minimum  residual  methods  has  been  observed  previously  [8];  however,  the  DPG  method 
distinguishes  itself  by  measuring  the  residual  of  the  operator  form  of  the  equation,  which  is  posed 
in  the  dual  space.  Soon  after,  Cohen,  Dahmen  and  Welper  introduced  an  alternative  saddle- 
point  formulation  of  such  a  minimum  residual  method  in  [2],  which  alluded  to  a  Variational 
Multiscale  (VMS)  perspective  [9,  10,  11]. 

The  goal  of  this  paper  is  three-fold.  The  first  is  to  present  a  method  which  borrows  concepts 
from  each  of  these  recent  works  and  to  derive  a  more  detailed  connection  between  the  Petrov- 
Galerkin,  stabilized,  and  VMS  perspectives  of  minimum-residual  methods.  The  second  is  to 
demonstrate  that  the  proposed  method  is  robust  in  both  convection  and  diffusion-dominated 
regimes  provided  the  dual  norm  used  to  measure  the  residual  is  chosen  intelligently.  The  last 
goal  is  to  show  that  the  method  is  stable  for  arbitrary  higher  order  and  adaptive  meshes,  easily 
implemented  using  existing  finite  element  codes  and  technologies,  and  extendable  to  three- 
dimensional  convection-diffusion  problems. 

2.  A  minimum-residual  method.  Our  starting  point  is  the  minimization  of  some  mea¬ 
sure  of  error  over  a  finite-dimensional  space.  We  begin  by  first  introducing  an  abstract  variational 
formulation 


f  Given  l  G  V*,  find  u  G  U  such  that 
{  b(u,v)=l(v ),  Vv  G  V, 


1We  note  that  one  cannot  recover  the  SUPG  formulation  beginning  with  the  strong  form  of  the  equations 
and  the  SUPG  test  functions;  the  interpretation  of  SUPG  as  a  Petrov-Galerkin  scheme  holds  only  locally,  on 
element  interiors. 
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where  &(•,•)  :  U  xV  — >  M  is  a  continuous  bilinear  form.  Throughout  the  paper,  we  assume 
that  the  trial  space  U  and  test  space  V  are  real  Hilbert  spaces,  and  denote  U*  and  V *  as  the 
respective  topological  dual  spaces.  Throughout  the  paper,  we  suppose  the  variational  problem 
(2.1)  to  be  well-posed  in  the  inf-sup  sense.  We  can  then  identify  a  unique  operator  B  :  U  V* 
such  that 


(Bu,v)y*xv  •=  b(u,v),  u  eU,v  eV 

with  (•,  • )y*xv  denoting  the  duality  pairing  between  R*  and  V,  to  obtain  the  operator  form  of 
the  our  variational  problem 

(2.2)  Bu  =  l  mV*. 


We  are  interested  in  minimizing  the  residual  over  the  discrete  approximating  subspace  Uh  C  U 


uh  =  argmin  J{uh)  := 

Uh^Uh 


\\\l-Buh\\2y 


1  \l(v)  -  HvjuV) |2 

2  rCl'\|ll)  My 


For  convenience  in  writing,  we  will  abuse  the  notation  supvGy  to  denote  supvGy^0}  f°r  the 
remainder  of  the  paper.  If  we  define  the  problem-dependent  energy  norm 


IMI E  ll^llv*  5 

then  we  can  equate  the  minimization  of  J(uh )  with  the  minimization  of  error  in  ||  'll#- 

The  first  order  optimality  condition  for  minimization  of  the  quadratic  functional  J(uh) 
requires  the  Gateaux  derivative  to  be  zero  in  all  directions  Su  G  Uh , 


(2.3)  (/  -  Buh,  BSu)v *  =  0,  G  U, 

which  is  nothing  more  than  the  least-squares  condition  enforcing  orthogonality  of  error  with 
respect  to  the  inner  product  on  V. 

The  difficulty  in  working  with  the  first-order  optimality  condition  (2.3)  is  that  the  inner 
product  (•,  -)y*  cannot  be  evaluated  explicitly.  However,  we  have  that 

(2.4)  (l  -  Buh,  B5u)v,  =  (Ry1^  -  Buh),RylB5u)v  , 

where  Ry  :  V  is  the  Riesz  map  mapping  elements  of  a  Hilbert  space  V  to  elements  of  the 

dual  R*  defined  by 


(RyV,  5v)v*  xV  ( v,5v)v  . 

Furthermore,  the  Riesz  operator  is  an  isometry,  such  that  J(uh)  =  \  ||Z  —  Buh\\y*  =  \  ||Ry1(/  —  Buh)\\y- 
Thus,  satisfaction  of  (2.4)  is  exactly  equivalent  to  satisfaction  of  the  original  optimality  condi¬ 
tions  (2.3). 

2.1.  Saddle  point  formulation.  We  can  translate  the  optimality  conditions  (2.4)  into  a 
more  standard  variational  problem.  First,  given  some  finite  dimensional  solution  Uh  G  C/^,  we 
define  the  error  representation  function  e  such  that 

(e,  v)v  =  l(v)  -  b(uh ,  v),  Vv  G  V. 

We  recognize  this  condition  as  defining  Rye  =  l  —  Buh\  inversion  of  the  Riesz  map  gives  us 

e  =  Ry1  ( l  -  Buh )  =  RyXB  (u  -  uh) , 

from  which  we  can  see  that  ||e||y  =  || u  —  Uh\\E- 

Secondly,  we  need  to  enforce  the  orthogonality  of  the  error  in  equation  (2.4),  which  becomes 

(e,  R^BSv^y  =  0,  \/Su  G  Uh 
— (^,  iBSujyxy*  —  b(5u,  e)  —  0. 
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Putting  these  two  conditions  returns  the  symmetric  saddle-point  problem  for  (e,Uh)  Uh 

(e,  v)v  +  b(uh,  v)  =  l(y),  Vu  G  V 
b(Su ,  e)  =  0,  WSu  G  Uh- 

We  note  that,  under  the  perspective  of  the  minimization  of  J(uh ),  this  system  can  be  viewed  as  a 
constrained  optimization  problem  for  e,  where  Uh  are  Lagrange  multipliers  enforcing  satisfaction 
of  the  variational  problem. 

At  the  moment,  V  remains  an  infinite-dimensional  space,  and  will  have  to  be  discretized  by 
some  Vh  C  V  in  order  to  yield  a  solution.  Under  this  infinite  dimensional  setting,  the  choice 
of  discretization  for  V  will  determine  how  effectively  the  solution  of  the  discrete  problem  will 
minimize  J(uh)  ~  inadequate  resolution  of  Vh  may  result  in  a  solution  which  is  far  from  the  true 
minimizer. 


(e,v)v 

- 

e 

b(5u ,  e) 

u 

Fig.  2.1.  Saddle  point  system  for  minimization  of  J(uy);  the  Schur  complement  eliminating  degrees  of 
freedom  for  e  returns  the  system  resulting  from  optimal  test  functions. 

This  formulation  was  first  introduced  by  Dahmen  et  al.  in  [2],  where  Vy  was  determined 
adaptively.  The  method  proceeded  in  two  steps  -  first,  Vy  was  set  equal  to  Uh  on  what  we 
will  define  as  the  coarse  mesh.  Next,  a  posteriori  error  estimators  for  e  (based  on  a  specific 
choice  of  inner  product  (’5  Oy)  were  used  to  drive  adaptivity  to  determine  a  fine-scale  mesh  on 
which  Vh  was  supported.  We  choose  a  simpler  representation  -  if  Uh  is  represented  by  piecewise 
polynomials  of  order  p,  we  set  Vy  to  be  piecewise  polynomials  of  order  p  +  A p.  We  note  that 
these  two  choices  of  discretization  for  V  are  not  mutually  exclusive,  and  that  novel  choices  for 
Vh  are  likely  the  key  to  yielding  computationally  effective  methods  under  this  framework. 

2.2.  DPG  optimal  test  function  framework.  The  DPG  method  was  introduced  by 
Demkowicz  and  Gopalakrishnan  in  2009,  and  rapidly  applied  to  a  series  of  problems  in  com¬ 
putational  mechanics  [1,  12,  13,  14,  15].  Though  the  DPG  method  shares  the  same  variational 
setting  as  Dahmen  et  al.  each  method  followed  very  different  approaches  in  implementation. 
Recall  the  orthogonality  of  error  (2.4);  by  defining  the  error  representation  function  e  and  un¬ 
doing  the  Riesz  map  acting  on  BSu ,  we  recovered  a  saddle  point  formulation.  However,  it  is 
possible  to  instead  undo  the  Riesz  map  acting  on  the  residual  instead,  to  recover 

[Ry1  ( l  -  Buh ) ,  R^BSi^y  =  0,  WSu  e  Uh 

=>  ((l  ~  Buh)  ^RVlB^U)V*xV  =  0 
=>■  l  (RyXB5u)  —  b  (uh^Ry1  B5u)  =  0 

If  we  define  vsu  :=  R^Bbu,  then  we  recover  the  original  variational  formulation  with  a  non¬ 
standard  test  space 


b(uh,vsu)  =  l(vsu ),  Vfe  G  Uh. 
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The  DPG  method  with  optimal  test  functions  was  constructed  under  this  perspective  of  mini¬ 
mization  of  J(uh ).  While  vsu  G  V  must  still  be  solved  for  approximately,  under  the  assumption 
of  discontinuous  test  functions,  this  solve  can  be  localized  to  a  single  element.2 

We  can  draw  yet  another  connection  to  the  DPG  method  by  considering  the  discrete  problem 
which  arises  from  the  saddle  point  formulation.  Our  discrete  saddle  point  system  is 


where,  for  Uj  G  £4,  and  Vi ,  Vj  G  W, 

-A-ij  •=  (vi ,  Uj ) y  ,  Bij  .=  b  (uj,Vi) , 

such  that  A  is  a  square  matrix  and  B  is  an  overdetermined  rectangular  matrix.  Expressing 
e  =  A-1  (/  —  Bu )  allows  us  to  statically  condense  out  the  degrees  of  freedom  e  to  form  the 
Schur  complement 

BTA1Bu  =  BTA1f. 

The  above  system  is  essentially  an  algebraic  least  squares  system,  but  preconditioned  on  the 
inside  with  the  positive-definite  operator  A,  which  is  precisely  the  discrete  system  that  arises 
under  the  optimal  test  function  framework  of  DPG.  We  can  see  this  further  by  noting  that 

(A~lB)T  Bu  =  (. A~lB)T  f 

is  equivalent  to 

/  dimVk  \  / dim  Vh  \ 

b  |  Ujn  ^  ^  Vik4>k  I  =  l  I  ^  ^  Vikcfrk  I  ,  i  —  1, . . . ,  dim(J7/l) 

where  vik  =  (A_1£?)^  are  the  degrees  of  freedom  associated  with  the  ith  approximated  optimal 
test  function. 

The  efficiency  of  the  DPG  method  lies  in  that,  under  the  assumption  of  a  localizable  norm 
and  optimal  test  functions,  A  is  block  diagonal  and  can  be  inverted  locally  [6].  Under  the 
ultra- weak  formulation  of  DPG,  one  can  show  that  this  corresponds  exactly  to  choosing  a  non- 
conforming  DG  space  for  V,  where  continuity  is  enforced  weakly  at  element  interfaces  by  re¬ 
quiring  that  the  jump  be  orthogonal  to  polynomials  of  order  p  [16]. 3 * 

2.3.  Variational  Multiscale  Framework.  The  above  saddle  point  system  can  also  be 
directly  derived  from  variational  multiscale  (VMS)  principles  [9,  10].  The  key  difference  between 
this  and  standard  multiscale  methods  is  the  way  in  which  the  fine  scales  influence  the  coarse; 
approximating  the  fine-scale  problem  is  as  difficult  as  approximating  the  original  equation, 
so  analytic  approximations  or  assumptions  must  typically  be  made.  However,  through  the 
proper  change  of  variables,  this  method  converts  the  fine-scale  problem  into  a  symmetric-positive 
definite  one,  allowing  for  a  well-behaved  subgrid  model  of  fine  scale  behavior. 

We  begin  again  with  the  variational  problem 

(Bu,v)y*xy  =  (l,v)y*xy  , 

where,  as  before,  B  :  U  V*  is  the  operator  form  of  b(u,v),  and  l  G  V*  is  the  load.  Standard 
VMS  techniques  then  decompose  U  into  a  coarse  scale  space  Uh  and  fine  scale  space  [T7,  such 

2In  practice,  if  Uh(K)  =  PP(K)  is  the  space  of  polynomials  of  order  p  over  the  element  K ,  then  the  enriched 
space  Vh  C  V  is  chosen  to  be  Vh(K)  =  Pp+Ap(i<r),  where  A p  >  n,  and  d  is  the  spatial  dimension. 

3Such  a  weakly  conforming  method  was  first  introduced  in  [17],  where  optimal  rates  of  convergence  were 

proved  under  the  assumption  that  the  DG  space  is  at  most  of  order  p  +  1  . 
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that  the  coarse  and  fine  scales  Uh  G  Uh  C  U  and  u!  G  U'  are  related  through  some  yet-unspecified 
orthogonality  condition.  Decomposing  iz  =  Uh  +  id  gives 

Uh  e  Uh,u'  e  u',v  e  V 

(BUfl,v)y,><y  +  ( BU  ,v)y,xy  =  (l ,  I})  y »  xy 

where  U  =  Uh®U'.  We  can  now  explicitly  choose  our  orthogonality  condition  defining  the  fine 
scales  and  rewrite  the  above  problem  as 

uh  eUh,u'  eU,v  eV 

{BUh,V)y,xy  +  ( Bll  ,V)y,xy  =  (l  ,  V)  y ,  xy 

(u',w)E  =  0,  \/weUh, 


where  ( u',w)E  is  defined  through 


( u  ,w) E  —  ( Tu  ,Tw)y 
T:U  ->  V,  T  :=  R^B 

where  (v,5v)v  is  an  inner  product  on  V  and  Ry  is  the  Riesz  map  associated  with  this  inner 
product.  DPG  literature  refer  to  T  as  the  trial-to-test  operator,  which,  given  a  basis  function 
( j>i  G  Uh ,  returns  the  optimal  test  function  for  [1]. 

If  we  define  our  fine-scale  using  the  above  orthogonality  condition,  we  have  that 

(Tu  , v)y  =  ( Ry  Bu  ,v)y  =  ( Bu  ,v)y*xy 
( u  =  ( Tu  ,  Ry  Bw^jy  =  (Tu  ,  Bw)yxy* 

and  can  now  rewrite  our  above  VMS  formulation  as 

uh  eUh,u'  eU,v  eV 

(BUh,v)y*  xy  +  (Tu  ,v)y  =  (l,v)y*  xy 

(Tu',  Bw)VxV *  =0,  \/w  G  Uh • 

The  saddle  point  formulation  of  Section  2.1  can  be  interpreted  as  solving  not  for  u! ,  but  for 
e  :=  TV.  The  discrete  approximation  then  follows  directly  from  setting  V  :=  Vh,  some  finite¬ 
dimensional  enriched  space,  which  impacts  the  method  only  through  the  approximation  of  Tu' 
and  the  construction  of  the  Riesz  map  Ry1 . 

To  summarize,  the  minimum-residual  framework  motivating  both  Dahmen’s  variational  sta¬ 
bilization  and  Demkowicz  and  Gopalakrishnan’s  DPG  method  can  be  characterized  as  a  VMS 
method  under  which 

•  Orthogonality  of  Uh  and  U'  is  defined  through  a  nonstandard  energy  inner  product  (i. 
e.  (u',uh)E  =  0),  and 

•  A  change  of  variables  (i.  e.  solving  for  Tu'  instead  of  u')  symmetrizes  the  fine-scale 
contribution  and  overall  system. 

3.  Choices  for  test  norms.  As  with  the  DPG  method,  the  primary  theoretical  difficulty 
is  in  defining  the  norm  |H|y  such  that  the  method  returns  good  results.  This  choice  of  norm 
defines  the  measure  of  the  residual  which  we  minimize  or  the  orthogonality  condition  defining  the 
separation  of  coarse  and  fine  scales,  depending  upon  whether  you  adopt  the  minimum-residual 
or  VMS  perspective. 

From  this  point  on,  we  will  focus  specifically  on  H 1  standard  Galerkin  formulation  for  the 
scalar  convection-diffusion  problem  on  domain  D 

V  ■  (Pu)  -  eAu  =  f, 
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which  represents  the  convection  of  some  quantity  u  subject  to  a  convection  field  (3  G  Rn  and 
diffusion  of  magnitude  e  G  M,  where  e  >  0.  We  impose  inflow  and  outflow  boundary  conditions 


u  I  Tout  -  U°ui 

where  the  following  boundaries  are  defined  as 

r in  :=  {x  G  dfls.  t.  (3  •  n(x)  <  0} 

Tout  :=  {%  G  dfls.  t.  (3  •  n(x)  >  0} 

where  n(x)  is  the  outward  normal  at  a  point  on  the  boundary  dfl.  If  it  is  non-empty,  we  may 
also  define  To  :=  {x  G  dQs.  t.  /3  •  n(x)  =  0}  as  the  boundary  with  outward  normal  orthogonal 
to  the  streamline. 

Under  m,vG  Hq(Q),  Our  weak  formulation  becomes 

(Vrt,  /3v  +  eVv) L2^  =  (f,v) L2^  • 

A  norm  based  on  the  symmetric  part  of  the  bilinear  form  was  explored  in  [2],  but  yielded  ap¬ 
proximations  which  degnerated  in  e.  We  present  here  two  different  methods  based  on  alternative 
choices  for  test  norms  on  V  for  the  convection-diffusion  equation  that  both  yield  good  results 
as  e  -A  0. 


3.1.  Weighted  streamline  diffusion  norm.  Our  first  method  is  motivated  by  previous 
results  for  the  DPG  method  in  [18],  which  achieved  good  results  through  a  modification  of  the 
test  norm.  We  set 


IMIv.tB  :=  ill  II •  V«|£a(n)  +  e  ||Vv||ia(n)  +  a  |Mlh(n) , 

where  L  =  |f2|3  is  the  length  measure  of  ft  in  d  spatial  dimensions,  \/3\  =  max^  ||/?||j2,  and  we 
have  introduced  the  scaling  ^  on  the  streamline  diffusion  term  for  unit  consistency4,  a  >  0  is  a 
user-specified  constant  with  units  of  inverse  time5,  and  w(x)  is  some  smoothly  varying  positive 
weighting  function.  Specifically,  we  require  \w(x)\  to  be  of  magnitude  0(e)  for  x  G  U n,  but  0(1) 
away  from  the  inflow  boundary.6 

A  second  motivation  behind  the  weighting  function  can  also  be  found  in  [19,  3,  5],  which 
refer  to  w(x)  as  a  cut-off  function.  A  theoretical  issue  with  the  convection- diffusion  problem  is 
that  | \u\ \Hi  and  \\u  —  Uh\\Hi  do  not  remain  bounded  as  e  — )>  0,  due  to  the  fact  that  u'  =  0(e_1) 
in  the  region  of  boundary  or  internal  layers.  These  cut-off  functions  were  introduced  in  order 
to  yield  e-independent  bounds  on  the  error  by  localizing  error  estimates  away  from  layers. 

We  can  demonstrate  the  need  for  the  weight /cut-off  function  w(x)  in  the  test  norm  -  Fig¬ 
ure  3.1  shows  the  result  of  two  simulations  on  a  uniform  16  element  piecewise-linear  mesh.  The 
error  representation  e  is  approximated  by  quadratics  over  the  same  mesh.  Without  the  weight, 
the  solution  degenerates  near  the  inflow;  a  similar  phenomenon  was  observed  in  [2,  1,  16]. 

Physically  speaking,  the  weight  can  be  interpreted  as  a  cutoff  function  for  the  adjoint  so¬ 
lution  [20],  which  displays  a  boundary  layer  at  the  inflow  due  to  the  fact  that  the  direction 
of  convection  is  reversed  for  the  adjoint  equation.  Mathematically  speaking,  it  can  also  be 
interpreted  as  restoring  a  sense  of  directionality  to  the  test  norm  (which  is  sign  of  /3). 


4In  all  problems  considered,  =  0(1). 

5In  the  case  where  a  =  0,  the  seminorm  remains  a  norm  due  to  equivalence  with  the  H 1  seminorm,  which 
is  norm  on  Hq.  Specifying  a  >  0  does  not  change  the  results  dramatically;  however,  it  does  appear  to  improve 
conditioning  of  A. 

6 Specific  restrictions  on  w  are  derived  in  [18]  in  context  of  the  DPG  method. 
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Fig.  3.1.  Comparisons  of  u  and  e  for  e  =  10  3,  f  =  1  and  a  =  0.  The  solution  when  w  =  1  is  in  red ,  and 
the  solution  for  when  w  =  x  +  e  is  given  in  blue.  The  exact  solution  is  given  in  black. 


3.2.  Streamline  norm  with  weak  boundary  conditions  on  e.  While  the  above  norm 
appears  to  give  improved  results,  the  issue  of  constructing  a  weighting  function  for  complex 
geometries  in  higher  dimensions  can  be  fairly  nontrivial.  One  possibility  would  be  to  solve  a 
Laplace’s  equation  for  w  on  the  domain,  enforcing  that  w(x)  =  e  at  inflow  boundaries;  however, 
the  behavior  of  such  a  computationally  determined  weight  has  not  been  examined,  and  incurs 
additional  complexity  that  is  avoided  by  other  methods. 

Under  this  saddle-point  framework,  however,  it  is  possible  to  avoid  the  use  of  a  weight 
altogether.  Currently,  we  seek  i?,  e  G  Hq(Q)]  however,  if  we  seek  instead 

V,ee  :=  {v  e  H\Q)s.  t.  v\Fout  =  0}  , 

we  modify  our  variational  formulation  to  be 


(Vu,Pv  +  eVv)L  2(n) 


du 

dnV 


(/^)l2(o)  • 


Under  this  modification  of  the  variational  formulation,  we  can  set  w  =  1  and  use  the  test  norm 


lit  :=  m  WP ' 


Vvllz,2(o)  +e||Vv||^2(n)  +a||v| 


2 

L2(0) 


without  experiencing  the  issues  displayed  for  e  G  Hq(Q)  in  Figure  3.1. 

Physically  speaking,  the  replacement  of  Hq(Q)  with  can  be  viewed  as  releasing  the 

strong  imposition  of  a  boundary  condition  at  the  inflow.  The  addition  of  the  term  e  fr  ^ v 
is  akin  to  a  penalty  term  that  weakly  enforces  a  homogeneous  boundary  condition  on  e  at  the 
inflow;  the  orthogonality  condition  for  e  under  this  additional  term  is 

(VSu,/3e  +  eVe)L2(n) e  f  =  0,  W5ueUh 

d  Tin 

If  Uh  is  taken  to  be  a  polynomial  space  of  order  p,  the  first  term  enforces  orthogonality  of  /?e+eVe 
to  polynomials  of  order  p  —  1  on  the  interior,  while  the  second  term  enforces  orthogonality  of  e  to 
polynomials  of  order  p—  1  on  the  subset  of  the  boundary  T[n.  If  e  is  represented  by  a  polynomial 
of  order  p—  1,  then  this  condition  reduces  to  strong  enforcement  of  e|r.  =  0  [17].  This  boundary 
condition  is  released  in  the  limit  as  e  0;  as  a  result,  for  small  e,  the  contribution  (e,  v)v  no 
longer  provides  an  extraneous  stabilization  at  the  inflow  for  under-resolved  meshes,  improving 
the  resulting  solution  u.  However,  for  large  e,  the  boundary  condition  activates  and  e  recovers 
homogeneous  boundary  conditions  over  the  entire  space  V.  The  choice  of  weight  w(x)  in  the 
previous  section  can  be  seen  as  similarly  removing  the  contribution  of  e  near  the  inflow  for  small 
e,  but  allowing  it  to  return  for  large  e. 

A  comparison  between  the  weighted  norm  IMIwu,  with  e  G  i^o(U)  and  unweighted  norm 
with  e  G  H^ut  in  one  dimension  is  given  in  Figure  3.2. 
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(a)  Solution  u 


(b)  Fine-scales  e 


Fig.  3.2.  Comparisons  of  u  and  e  for  e  =  10  1 
e  E  H[ut(Q)  is  in  red ,  and  the  solution  for  when  w  =  x- 


,  f  =  1  and  a  =  0.  The  solution  when  w  =  1  and 
-  e  is  given  in  blue.  The  exact  solution  is  given  in  black. 


Remark  1.  We  remark  that  this  inclusion  of  e  fr  is  a  variational  crime  under  the  assump¬ 

tion  that  u  G  Hq(D);  the  proper  way  to  include  such  a  term  in  a  conforming  manner  (without 
introducing  additional  regularity  assumptions  on  u)  would  be  to  introduce  as  an  additional 
unknown  on  the  boundary,  similar  to  the  work  done  using  DPG  under  the  ultra-weak  formulation 
[1,  21].  We  note  that  Broersen  and  Stevenson  have  also  achieved  good  results  by  making  use  of 
the  space  H[ut(D)  in  their  application  of  the  DPG  method  to  the  mixed  mild-weak  formulation 
for  convection- diffusion  [22].  The  mild-weak  formulation  effectively  enforces  the  orthogonality 
of  e  to  polynomials  of  lower  order  on  rin  in  a  conforming  manner,  by  introducing  an  additional 
flux  unknown  representing  e|^. 

For  many  flow  problems,  activity  near  the  inflow  is  typically  very  smooth  and  can  thus  be 
represented  using  a  minimal  number  of  elements  under  an  effective  adaptive  mesh  refinement 
scheme,  implying  that  the  cost  of  introducing  such  additional  unknowns  on  T in  can  be  minimized. 
However,  since  the  codebase  we  currently  use  does  not  feature  unknowns  with  support  only  on 
element  edges/faces,  we  continue  to  operate  under  this  variational  crime,  in  hopes  that  we  can 
implement  a  conforming  version  in  the  near  future. 

Remark  2.  Various  other  methods  of  implementing  a  weak  boundary  condition  on  e  at  the 
inflow  T in  were  implemented ;  for  example,  the  inner  product  ( e,v)v  was  modified  by  adding  the 
a  weak  penalization  term 


eC(h,p 


ev 


to  enforce  the  boundary  condition  at  e|r.n,  where  C(h,p+Ap)  =  p2/h,  which  is  a  common  choice 
of  penalization  parameter  in  high  order  finite  element  literature  [23,  24].  However,  the  inclusion 
of  the  nonconforming  term  e  fr  Qffv  in  the  bilinear  form  still  resulted  in  the  best  solutions  for 
all  our  numerical  tests. 

Thus,  the  full  formal  characterization  of  the  second  method  can  be  given  as  follows  -  solve 
for  u  euh  <Z  Hq(Q)  n  Hz+6(Q)  (where  S  is  an  positive  real  number),  e  G  Vh  C  H[ut(Q)  such 
that 


f  du 

(e,v)v  +  (Vu,/3v  +  eVv)L 2(n)  -  e  /  —v  =  (/,  v)L2{n),  \/v  G  Vh 

^  Tin 


dSu 


(\75u,[3e  +  eVe)L2^  -  ej  ~^e  =  0,  VSueUh 


where 


(e,  v)v  •  Ve,  /3  •  Vu) L2^  +  e  (e,  Vu) l2 ^  • 


We  will  utilize  this  second  version  of  the  method  for  all  following  numerical  examples. 


Remark  3.  We  note  that ,  when  the  Peclet  number  Pe  =  >>  1,  the  boundary  term  e  fr 

effectively  disappears ,  and  numerically  speaking ,  we  can  use  e, v  G  H^ut(Q)  without  including 
the  additional  boundary  term  in  the  variational  formulation. 

4.  Numerical  examples.  We  present  several  numerical  examples,  demonstrating  rates  of 
convergence  and  performance  of  the  method  on  several  benchmark  problems  in  2D.  All  numerical 
examples  were  implemented  in  the  Python-based  FEniCS  codebase  [25]  to  demonstrate  that, 
unlike  previous  methods  based  upon  the  minimization  of  dual  residuals,  this  method  is  easily 
implemented  in  most  standard  high-order  finite  element  codes. 

A  remaining  choice  to  be  specified  is  the  discretization  of  VA  if  Uh  is  given  to  be  polynomials 
of  order  p,  the  enriched  space  Vh  under  which  the  fine-scales  are  approximated  is  set  to  be 
polynomials  of  order  p+ Ap  -  in  other  words,  we  use  a  simple  polynomial  enrichment  scheme  over 
a  fixed  mesh  to  approximate  the  fine-scale  error  representation  function  e.  In  all  experiments, 
A p  =  1  has  been  found  to  be  sufficient,  and  larger  A p  has  not  been  found  to  have  a  significant 
effect  on  the  behavior  of  the  method.  This  is  in  contrast  to  the  original  DPG  method,  which 
also  adopts  a  polynomial  enriched  space  (where  enrichment  is  done  locally  over  each  element) 
-  results  by  Gopalakrishnan  and  Qiu  in  [26]  suggest  that  only  Ap  >  d ,  where  d  is  spatial 
dimension,  is  guaranteed  to  yield  optimal  convergence  rates  for  Laplace’s  equation. 

We  report  only  L 2  errors  for  each  example,  due  to  the  fact  that  H 1  error  is  not  a  well- 
defined  quantity  as  e  — 0.  As  a  consequence,  during  adaptive  refinement,  the  error  does  not 
always  monotonically  decrease  or  behave  as  expected,  since  adaptive  refinement  is  driven  by  the 
the  quantity  ||  u  —  Uh  We  =  IMIv.  which  we  expect  to  be  equivalent  (possibly  with  equivalence 
constants  depending  on  e)  to  the  H 1  norm.  We  do  expect  || u  —  Uh\\E  to  decrease  monotonically 
for  all  e  (though  not  uniformly  in  e)  due  to  the  minimum-residual  nature  of  our  method. 

4.1.  Erikkson- Johnson  model  problem.  We  adopt  a  problem  first  proposed  by  Eriks¬ 
son  and  Johnson  in  [27]  and  later  used  in  [18,  20]  to  determine  the  robustness  of  DPG  with 
respect  to  the  diffusion  parameter  e.  For  the  choice  of  Q  =  (0,  l)2,  /  =  0,  and  (3  =  (1,0)T,  the 
convection  diffusion  equation  reduces  to 

du  ( d2u  d2u\ 

dx  6  \  dx 2  dy2  J 

which  has  an  exact  solution  by  separation  of  variables,  allowing  us  to  analyze  convergence  of 
DPG  for  a  wide  range  of  e.  For  boundary  conditions,  we  impose 

u  -  Uq:  x  =  0, 
u  =  0,  x  =  1,  y  =  0, 1 

In  this  case,  our  exact  solution  is  the  series 

eri,n(x-l)  _  er2,n(x- 1) 

u(x,  j)aV  Cn - — - — - sm(niry) 

f  ^  g  '  1  ,n  —  g  '  2,n 

n=0 


where 


rl,2,n  = 


i  ±  Vi  +  4An 

2e 


\  2  2 
\n  =  n  7r  e. 

We  take  for  the  exact  solution  the  first  mode,  such  that  Co  =  1,  and  Cn  =  0  for  n  1. 

en(x-l)  _  er2(x- 1) 


u(x,y)  = 
r  1,2 


e~ri  -  e~r2 
1  zb  y/l  +  4e27T2 


sin(7n/) 


2e 
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The  problem  is  driven  solely  by  inflow  boundary  conditions,  and  develops  a  boundary  layer 
of  width  e  at  the  outflow  x  =  1.  The  resulting  solution  is  shown  in  Figure  4.1  for  e  =  10-2, 
along  with  convergence  rates  of  L2  (Q)  error  under  uniform  mesh  refinement  for  both  p  =  1  and 
p  —  2.  Optimal  rates  of  convergence  are  not  expected  due  to  the  pre- asymptotic  nature  of  the 


(a)  Exact  solution  for  u  (b)  Convergence  rates  for  e  =  10  2 

Fig.  4.1.  Erikkson- Johnson  model  problem  for  e  =  10— 2 ,  along  with  convergence  rates  under  uniform 
refinement. 


mesh  relative  to  the  solution  -  typically,  h  must  be  O(e)  in  order  to  observe  optimal  rates  of 
convergence  [28]. 

For  large  diffusion  -  e  =  1.0,  we  do  observe  optimal  p+  1  rates  of  convergence  in  Figure  4.2. 
Likewise,  when  e  «  h,  we  observe  a  sub-optimal  rate  of  |,  the  same  as  the  rate  suggested  by 
theory  and  observed  numerically  in  [22]  due  to  the  strong  boundary  layer  present  under  small 
e. 


Fig.  4.2.  Erikkson- Johns  on  model  problem  for  e  =  1  and  e  =  10  4;  along  with  convergence  rates  under 
uniform  refinement. 
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5.  Adaptive  mesh  refinement.  We  experimented  also  with  an  adaptive  refinement  scheme. 
Since  IMIv  is  localizable7,  we  can  evaluate  it  element-wise  to  get  element  error  indicators 
ck  :=  llelly( k) •  We  implemented  a  bulk-chasing  refinement  strategy,  where,  given  some  fac¬ 
tor  6  G  [0, 1],  we  refine  the  top  \0N]  elements  with  the  largest  error  indicators  e^.  A  greedy 
refinement  scheme  was  also  implemented,  where  all  elements  with  >  9  max^  eE  are  refined; 
however,  this  refinement  scheme  tended  to  place  refinements  more  conservatively,  requiring 
many  more  refinement  iterations  to  achieve  a  qualitatively  good  resolution.  In  all  following 
experiments,  6  is  set  to  be  .25. 

5.1.  Erikkson- Johnson  model  problem.  We  continue  to  verify  our  method  using  the 
Erikkson- Johnson  with  e  =  10-3.  The  solution  u  and  the  convergence/rates  of  L 2  (Cl)  error 
under  adaptive  mesh  refinement  for  p  =  2  are  given  in  Figure  5.1.  The  energy  error  \\u  —  Uh\\E 
is  also  displayed  for  comparison. 


Fig.  5.1.  Erikkson- Johns  on  problem  with  e  =  10  3  under  16  refinements. 

As  a  result  of  the  saddle  point  formulation,  we  can  also  examine  the  fine-scale  error  rep¬ 
resentation  function  e.  Figure  5.2  displays  a  comparison  between  the  error  representation 
function  after  the  first  iteration  of  adaptive  mesh  refinement  and  the  16th  iteration.  Recall 
that  the  main  contribution  to  the  energy  error  is  the  streamline  derivative  of  e;  \\u  —  Uh\\2E  = 
IMIv  =  ll/3-Ve||F([2)  +  e||Ve||F([2)-  Thus,  variation  in  the  streamline  direction  is  picked  up  by 

|| (3  •  Ve||^2(o)  and  is  the  primary  contribution  to  the  total  error.  After  sufficient  refinement  and 
resolution  of  the  outflow  layer,  the  error  representation  function  becomes  smooth  in  the  region 
of  the  boundary  layer,  and  refinements  will  be  placed  according  to  the  viscous  error  e  ||Ve|| 2L2^y 

One  open  question  concerning  this  method  is  the  energy  norm  \\u\\E  which  is  induced  by 
our  choice  of  bilinear  form  and  inner  product  on  our  choice  of  test  space  V  =  H^ut.  Figure  5.3 
displays  both  L 2  and  energy  errors  over  a  range  of  e.  We  observe  that  the  energy  error  appears 
to  be,  like  the  L 2  (Cl)  error  (but  unlike  H1^)  or  error  in  the  streamline  diffusion  norm),  inde¬ 
pendent  of  e  -  convergence  does  not  appear  to  be  in  a  norm  that  is  equivalent  to  the  H1  norm 
with  constants  independent  of  e. 

Numerical  experiments  indicate  that  the  measure  of  energy  error  is  fairly  insensitive  to 
the  degree  of  enrichment  of  W,  so  we  expect  that  only  a  moderately  finer  resolution  of  Vh 


7 A  localizable  norm  IMIy^)  can  be  written  in  the  form 

II  II 2  _  II  II 2 

IMIv(Q)  —  IMI V(K)  ’ 

Kenh 

where  |M|y(x)  a  norm  over  each  element  K  in  the  mesh 
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Fig.  5.2.  Comparison  of  the  fine  scale  error  representation  e  under  the  Erikks on- Johnson  problem  with 
e  =  10-3  for  the  mesh  at  the  first  and  1 6th  refinement  steps. 


with  respect  to  Uh  is  sufficient  to  capture  fine-scale  effects.  More  investigation  is  necessary  to 
determine  the  nature  of  the  energy  norm  induced  by  this  method. 


Fig.  5.3.  Comparison  of  L2  (Q)  and  energy  error  convergences  under  adaptive  refinement  for  the  Erikkson- 
Johnson  problem  for  e  =  10-3, 10  4 , 10-5, 10-6. 


5.2.  Discontinuous  forcing  and  advection  skew  to  mesh.  We  consider  next  a  slight 
variant  of  a  model  problem  introduced  by  Hughes  and  Sangalli  in  [11].  The  domain  Cl  is  again 
taken  to  be  the  unit  square,  j3  is  taken  to  b e  /3  =  [.5, 1],  and  e  =  10-3.  Homogeneous  Dirich- 
let  boundary  conditions  are  taken  over  the  entire  boundary,  and  the  problem  is  driven  by  a 
discontinuous  forcing  term  f(x,y),  where 


fj,y) 


l\iy>2x 
0  otherwise. 
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Under  such  an  example,  a  boundary  layer  develops  at  the  outflow  boundaries  x,y  =  1,  and  an 
internal  layer  develops  due  to  the  fact  that  the  discontinuity  in  /  is  parallel  to  the  streamline. 
Figure  5.4  displays  the  solution  and  mesh  resulting  from  16  automatic  mesh  refinement  itera¬ 
tions.  Figure  5.5  displays  a  zoom  of  the  solution  and  convergence  history  of  the  energy  error. 


Fig.  5.4.  Model  problem  with  e  =  10  3  and  discontinuous  forcing  f. 


We  can  similarly  examine  advection  skew  to  the  mesh,  another  common  benchmark  problem 
for  convection-diffusion.  We  adopt  a  modification  of  the  standard  problem  parameters  -  here, 
e  =  10-7,  (3  =  [.5, 1],  and  f  —  0.  Boundary  conditions  are  set  such  that 


u  — 


1  if  y  =  0,  x  <  .5 
1  —  y  if  x  =  0 
0  otherwise. 


Mathematically  speaking,  the  considered  problem  is  ill-posed  for  the  standard  Galerkin  method 
in  the  continuous  setting  -  the  boundary  condition  implies  a  discontinuity  in  u ,  such  that  the 
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solution  cannot  lie  in  H1 .  The  boundary  condition  is  then  convected  in  a  manner  that  is  skew 
to  the  mesh  and  forms  a  boundary  layer  on  part  of  the  outflow. 


(a)  Solution  u 


(b)  u 


(c)  Error  representation  e 


Fig.  5.6.  Advection  skew  to  mesh  on  a  uniform  quadratic  32  x  32  grid  of  triangles.  The  fine  scale  error 
representation  function  e  is  also  shown;  note  that  the  multiscale  contribution  of  e  is  dependent  on  the  form  of 
(e,v)v,  which  involves  streamline  derivatives. 


Fig.  5.7.  Advection  skew  to  mesh  after  iterations  18  of  adaptive  mesh  refinement. 


We  note  that  the  standard  boundary  condition  for  advection  skew  to  the  mesh  is 


u  = 


lifx  =  0oiy  =  0  and  x  <  .5 
0  otherwise, 


which  produces  a  discontinuity  at  the  upper  left  hand  corner  as  well.  However,  under  automatic 
mesh  refinement,  extraneous  refinements  are  placed  at  this  corner’s  discontinuity.  It  is  possible 
to  avoid  such  an  issue  either  by  tailoring  a  mesh  refinement  scheme  to  avoid  the  corner  or  by 
adopting  a  regularized  boundary  condition  by  introducing  a  mesh-dependent  “ramped”  bound¬ 
ary  condition  which  approximates  the  discontinuity.  This  was  done  for  for  the  lid-driven  cavity 
flow  problem  under  Stokes  flow  in  [15]. 

We  avoid  both  the  above  “fixes”  in  order  to  more  cleanly  demonstrate  the  behavior  of  auto¬ 
matic  adaptive  mesh  refinement;  however,  we  note  that  under  the  uniform  meshes  in  Figure  5.9, 
or  under  adaptive  refinement  schemes  that  avoids  the  corner  singularity,  the  method  appears  to 
deliver  similar  results  to  SUPG. 
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Fig.  5.8.  Energy  error  for  advection  skew  to  mesh  over  18  iterations  of  adaptive  mesh  refinement. 


(a)  p  =  1  (b)  V  =  2 


Fig.  5.9.  Solution  u  for  advection  skew  to  mesh  with  the  upper-left  hand  corner  discontinuity  under  both 
linear  and  quadratic  uniform  meshes. 


5.3.  A  3D  model  problem.  We  conclude  by  applying  the  method  to  a  simple  model 
problem  in  3  space  dimensions.  We  solve  on  the  unit  cube  =  [0,  l]3  with  (3  =  [1,1,1],  /  =  i 
and  u  =  0  on  <9fh  A  boundary  layer  develops  on  the  outflow  boundaries  x,y,z  =  1.  Figure  5.10 
displays  slices  of  the  domain  along  the  streamline  direction  to  illustrate  the  formation  of  such  a 
boundary  layer. 

We  examined  the  behavior  of  the  method  with  respect  to  a  small  e,  in  this  case  e  =  10-6. 
8  adaptive  refinements  were  done  beginning  with  a  uniform  linear  mesh  of  8  elements  per  side; 
Figure  5.11  shows  the  resulting  solution  and  mesh  along  a  cut  orthogonal  to  the  streamline 
direction.  All  experiments  indicate  that  the  method  appears  to  behave  qualitatively  the  same 
as  in  2  dimensions,  modulo  increased  computational  cost. 

6.  Conclusions  and  future  direction.  We  have  presented  a  higher  order  adaptive  H1- 
conforming  method  for  convection-diffusion  problems  for  very  small  e,  and  have  made  connec¬ 
tions  to  both  the  DPG  method  as  well  as  the  method  of  Cohen,  Dahmen  and  Welper.  The 
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Fig.  5.10.  3D  Model  problem  with  e  =  5e  —  2. 


Fig.  5.11.  3D  Model  problem  with  e  =  10— 6  after  8  refinements.  The  cube  is  sliced  along  a  plane  with 
normal  direction  orthogonal  to  the  streamline  in  order  to  illustrate  the  clustering  of  refinements  in  the  interior 
of  the  domain. 


method  has  also  been  shown  to  be  derivable  in  the  Variational  Multiscale  framework  as  well, 
and  is  distinguished  from  traditional  VMS  methods  in  its  definition  of  the  fine-scale  space  and 
problem.  Unlike  standard  stabilized  H 1  methods,  there  are  no  stabilization  parameters.  How¬ 
ever,  unlike  the  DPG  methods  of  both  Demkowicz  and  Gopalakrishnan,  as  well  as  the  method 
of  Broersen  and  Stevenson,  this  method  is  easily  implementable  within  current  finite  element 
codes  in  arbitrary  spatial  dimensions. 

We  note  that,  while  the  method  displays  similarities  to  SUPG,  we  do  not  observe  nodal 
exactness  in  ID.  The  two  methods  differ  in  that  we  do  not  work  directly  with  coarse  scales 
defined  by  the  H 1  inner  product;  unless  we  can  induce  an  energy  inner  product  such  that 
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(■ Uh,u')E  =  (uh,uf)Hi,  then  we  should  not  expect  nodal  exactness  in  ID.  However,  in  higher 
space  dimensions,  the  behavior  of  the  method  is  qualitatively  very  similar  to  SUPG. 

We  hope  to  explore  in  future  research  both  additional  improvements  on  the  method  and 
the  application  of  this  method  to  other  problem  frameworks  and  settings,  a  few  of  which  are 
described  in  the  following  sections. 

6.1.  hpr- refinement.  Since  this  method  is  derived  from  the  same  variational  principles 
as  DPG,  it  should  likewise  be  automatically  inf-sup  stable  for  arbitrary  hp  meshes  as  well,  hp- 
refinement  is  of  interest  due  to  the  fact  that,  under  proper  choices  between  h  and  p  refinement, 
exponential  convergence  rates  are  possible,  and  have  been  observed  for  several  model  problems 
[29].  Intelligent  hp-adaptive  strategies  could  also  offer  more  mathematically  motivated  alterna¬ 
tives  to  CFD  approaches  such  as  flux  limiters,  which  effectively  limit  the  order  of  approximation 
in  areas  of  singularities  and  high  gradients. 

Additionally,  the  method  automatically  extends  to  incorporate  anisotropic  elements  and 
r-refinement  as  well,  similar  to  ALE  methods.  The  combination  of  these  two  methods  has 
been  shown  to  produce  very  fine-resolution  results  at  a  competitive  computational  cost  for  flow 
problems  with  very  small  viscosity  [30]. 

6.2.  Computational  feasibility.  The  main  obstacle  to  making  this  method  computa¬ 
tionally  competitive  is  the  doubling  of  the  number  of  fully  coupled  degrees  of  freedom  in  solving 
the  problem  in  saddle  point  system  formulation.  While  a  larger  number  of  degrees  of  freedom 
is  often  cited  as  a  detriment,  it  is  often  accompanied  by  an  advantageous  underlying  structure 
-  for  example,  discontinuous  Galerkin  methods  are  also  often  criticized  for  their  “explosion  of 
degrees  of  freedom”  compared  to  continuous  Galerkin  methods  [31],  but  their  local  nature  aids 
significantly  in  parallelization,  adaptivity,  and  explicit  time  integration,  as  opposed  to  methods 
based  on  continuous  discretizations. 

The  cost  of  explicitly  discretizing  the  fine  scale  error  representation  eG  roughly  amounts 
to  the  cost  of  discretizing  an  additional  coupled  PDE  under  a  higher  resolution.  The  idea  of 
solving  an  additional  PDE  for  the  purposes  of  stabilization  is  not  unheard  of  [32];  however, 
the  additional  discretization  cost  must  be  minimized.  We  propose  several  preliminary  ideas  for 
doing  so: 

•  Recall  in  Section  2.1  that  the  saddle  point  system  is  reducible  to 

BTA1Bu  =  BTA1f . 

where  /  and  B  are  the  load  and  overdetermined  stiffness  matrix  under  W,  and  A  is 
the  Gram  matrix/discrete  Riesz  map  associated  with  the  inner  product  (v,Sv)v.  One 
approach  is  to  explore  possibilities  for  rapid  inversion  of  A.  For  example,  it  is  possible 
to  choose  Vh  to  simply  be  the  space  of  piecewise  linears  enriched  with  bubbles  -  higher 
order  degrees  of  freedom  could  be  condensed  out,  resulting  in  a  system  for  e  with  the 
same  number  of  unknowns  as  u.  The  resulting  matrix-vector  products  required  to  form 
the  Schur  complement  could  be  parallelized  as  well. 

One  additional  possibility  would  be  the  exploration  of  iterative  methods  for  the  Schur 
complement  -  the  matrix  is  positive  definite,  and  assembling  a  matrix-vector  product 
could  be  done  efficiently  under  a  good  preconditioner  for  the  symmetric  positive-definite 
matrix  A-1.  Preliminary  results  demonstrate  effective  speedup  for  uniform  meshes  of 
piecewise  linears  in  2D  using  standard  preconditioners  (ILU,  block  Jacobi,  AMG),  and 
suggest  that  a  discretization  of  Vh  using  a  h-refined  fine  mesh  could  be  possible;  however, 
more  work  would  be  required  to  develop  preconditioners  for  high  order  adaptive  meshes. 

•  Another  possibility  would  be  to  solve  the  saddle  point  problem  using  an  Uzawa  method 
[33],  which  would  decouple  the  saddle  point  system  into  two  problems,  each  involving 
only  degrees  of  freedom  for  e  or  for  u.  The  solution  of  these  two  problems  feeds  an 
iterative  procedure  which  can  converge  to  the  solution  of  the  saddle  point  problem. 
Again,  effective  use  of  an  Uzawa  algorithm  would  require  a  good  preconditioner  for  A 
under  high-order  refined  meshes. 
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•  We  hope  to  explore  alternative  discretizations  for  both  Uh  and  Vh  in  order  to  minimize 
degrees  of  freedom  or  solution  for  in  the  overall  system.  A  few  examples  of  nonstan¬ 
dard  discretizations  include  higher  order  continuity  basis  functions  (splines  and  NURBS 
[34]),  and  discontinuous  functions  (DG).  Preliminary  experiments  have  yielded  promis¬ 
ing  results  in  ID  under  the  combination  of  higher-continuity  and  Co  bases  for  Uh  and 
W,  respectively.  We  hope  to  further  explore  intelligent  mixing  of  discretizations  for 
trial  and  test  spaces  in  future  experiments. 

•  The  DPG  method  achieves  computational  efficiency  by  block-diagonalizing  A  using 
discontinuous  test  functions.  The  same  method  could  be  applied  to  the  H 1  setting  for 
convection-diffusion  problems  at  the  cost  of  introducing  inter-element  fluxes  represent¬ 
ing  the  trace  on  the  boundary  of  an  element  resulting  from  integration  by  parts  of  the 
viscous  term.  An  alternative  along  the  same  lines  would  be  to  consider  block  diagonal- 
ization  of  A  under  a  patch-based  decomposition,  such  as  FETI  or  one  of  its  variants 
[35,  36]. 
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